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Abstract 

In this note, we show that there exist solutions of the Muskat problem that shift stability 
regimes: they start unstable, then become stable, and Hnally return to the unstable regime. We 
also exhibit numerical evidence of solutions with medium-sized L°° norm of the derivative of the 
initial condition that develop a turning singularity. 
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1 Introduction 

In this paper, we study two incompressible fluids with the same viscosity but different densities, 
and p~, evolving in a two dimensional porous medium with constant permeability k. The velocity v 
is determined by Darcy’s law 

= ( 1 . 1 ) 

where p is the pressure, p > 0 viscosity, and g > 0 gravitational acceleration. In addition, v is 
incompressible: 

V-t; = 0. (1.2) 

By rescaling properly, we can assume p = g = 1. The fluids also satisfy the conservation of mass 
equation 

dtp + v-yp = 0. (1.3) 

This is known as the Muskat problem m- We denote by D+ the region occupied by the fluid 
with density p'^ and by the region occupied by the fluid with density p~ ^ p'^. The point (0, oo) 
belongs to 12+, whereas the point (0, —oo) belongs to fl~. All quantities with superindex ± will refer 
to D+ respectively. The interface between both fluids at any time f is a curve z{a,t). We will work 
in the setting of flat at infinity interfaces, although the results can be extended to the horizontally 
periodic case. 

A quantity that will play a major role in this paper is the Rayleigh-Taylor condition, which is 
defined as 


RT{a, t) = - [Vp (z(a, t)) - Vp+(z(a, f))] • d^z{a, t), 

where we use the convention {u,v)^ = {—v,u). If RT{a,t) > 0 for all a S K., we will say that the 
curve is in the Rayleigh-Taylor stable regime at time t, and if RT{a,t) < 0 for some a G R, we will 
say that the curve is in the Rayleigh-Taylor unstable regime. 

One can rewrite the system (ITT])-(f01) in terms of the curve z(a,t), obtaining 

dtz{a,t) = ^ ^ ^ — P^7p{daz{a,t)-dpz{(3,t))dfl. (1.4) 

27r J^\z{a,t) - z[/3,t)\^ 
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A simple calculation of the Rayleigh-Taylor condition in terms of z{a,t) yields 


RT{a,t) = g{p - p~'')daZi{a,t). 


When the interface is a graph, parametrized as z{a,t) = (a,/(a, t)), equation (11.41) becomes 

o+ r 


dtf{x,t) = 


P - P 
27r 


x-y 


{x-yf+ {f{x,t)- 


{da,f{x,t) - dyf{y,t))dy 


(1.5) 


and the Rayleigh-Taylor condition simplifies to 

RT{a,t) = g{p~ - p+). 


The curve is now in the RT stable regime whenever p'^ < p~, i.e., the denser fluid is at the bottom. 

The Muskat problem has been studied in many works. A proof of local existence of classical 
solutions in the Rayleigh-Taylor stable regime and ill-posedness in the unstable regime appears in [9] . 
A maximum principle for ||/(■,t)||L“ can be found in [10) . Moreover, the authors showed in m that 
if ||i9x/o||l“ < 1; then \\dxf{-,t)\\L°° < ||9a;/o||L~ for all t > 0. Further work has shown existence of 
finite time turning (i.e., the curve ceases to be a graph in finite time and the Rayleigh-Taylor condition 
changes sign) [^. The precise result is the following: 


Theorem 1.1 There exists a nonempty open set of initial data in with Rayleigh-Taylor strietly 
positive (i.e., RT{a,0) > 0 for all a € sueh that the solutions of (11.41) have RT{a,t) < 0 for some 
finite time t > 0 and all a in some nonempty open interval. 

After this shift of regime, the curve may lose regularity. This was proved in [3]. More general 
models which take into account finite depth or non-constant permeability that also exhibit turning 
were studied in mis], where the estimates in the latter one were carried out by rigorous computer- 
assisted techniques, as opposed to the traditional pencil and paper ones from the former. All these 
results are local in time, and the techniques employed to prove them rely on a local analysis of 
the equations, therefore it is not possible to conclude global properties of the solutions from them. 
Concerning global existence, the first proof for small initial data was carried out in in the case 
where the fluids have different viscosities and the same densities (see also [5] for the setting of this 
paper: different densities and the same viscosities). Global existence for medium-sized initial data 
was established in Ii[7]. In the case where surface tension is taken into account, global existence was 
shown in [mill]. Global existence for the confined case was treated in [14]. Recently, in |^, a new 
framework was used to study global existence. The following theorem was proved in jSj: 


Theorem 1.2 Suppose that ||/o||l°° < oo and ||9x/o||l°° < l.Then there exists a global in time weak 
solution of CH) that satisfies 

f e C'([0,T] X R) nL°“([0,r];IF^’°°(IR)). 

In particular, f is Lipschitz continuous. 

The question of whether the constant 1 (which is independent of the physical parameters of the 
system) in the L°° bound of the derivative of the initial data is sharp or not is still open. The proof 
of Theorem 11.11 constructs initial conditions that have parameters that need to be taken big enough 
to ensure turning. Taken both these two facts into account, our initial motivation was to quantify 
and bridge the gap, finding or approximating the threshold C for the norm of dxfo such that 
initial data for which ||9 x/o||l°“ < C exist globally, and for each e > 0 there exist solutions with 
||5 x/o||l°“ = C -I- e which develop turning singularities. 

This note is organized as follows: in Section [3] we describe the numerical experiments that became 
the motivation for our main theorem, which we prove in Section |3| Finally, in Section |3| we discuss 
open problems and future work. 
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2 Numerical results 


In this section we describe numerical experiments which led to our main result, Theorem 13.31 below. 

Our motivation was to find initial data, as small as possible, which develop turning in finite time. 
For simplicity, we performed simulations in the horizontally periodic scenario. The evolution equation 
for the interface reads in that case: 


dtz{a,t) = 


47r 


iz{a,t) - z{/3,t)) 


sm(zi{a,t) - zi{/3,t)) 


cosh(z 2 (Q;,t) - Z2{l3,t)) - cos{zi{a,t) - zi(/3,t)) 


dp. ( 2 . 1 ) 


Picking different initial conditions and evolving them until they either develop turning or flatten out 
seems like looking for a needle in a haystack. Instead, we took data which we were sure that would 
turn (the interface has a vertical tangent at a point and the velocity is pointing in the right direction) 
and run the equation backwards in time for a short time to find our desired initial condition. 

At the linear level and backwards in time, equation (EH) behaves like a backwards heat equation, 
which is ill posed if the lighter fluid is on top of the denser one. To overcome this difficulty, we use 
the following heuristic: we perform very small steps backwards in time and at every step, we smooth 
the function by eliminating all frequencies whose components are below a given threshold. The reason 
behind this heuristic is that the family of solutions which turn over is an open set, and therefore 
small perturbations (the regularized versions) of the solution should remain in it. We remark that the 
backwards evolution is not done with the purpose of finding a numerical solution of the equation, but 
only to gain intuition about what initial condition to choose. Once we find a suitable candidate, we 
check its validity by running the equation forward with the candidate. 

The smoothing threshold was set to 10“®. The time integration was done using a Runge-Kutta 45 
scheme. The derivatives were calculated using a spectral method with a cutoff filter given by 


p{k) = exp 


-10 




In order to evaluate the singular integrals, we used an alternating quadrature rule [Dill]: 


dtz{ai,t) ss 2h- 


47r 




sin(2;i(ai,t) - zi(aj,t)) 


j—i odd 


cosh(z 2 (aj,t) - Z 2 {aj,t)) - cos{zi{ai,t) - zi{aj,t))' 


where h = ^ and N is the total number of nodes. We chose the condition at time t = 0 to be 


zi{a, 0) = a — sin(a). 


. 3sin(a) + 8sin(2a) + 3sin(3a) 

Z2(a,0) = --- 


and N = 2048 nodes (see Figured])- 

The smoothing was done after every At = 4 • 10“® and the simulation ran until t/ = —4.92 ■ 10“^. 
The obtained evolution of z{a,t) is depicted in Figure El 

Finally, we computed the evolution of equation (EH forward in time taking as initial condition 
z{a, —4.92-10“^), which has two vertical tangents at approximately (±3.795-10“^, ±1.268). The initial 
data is now in the stable regime and the equation is well-posed. We can use the same integration 
scheme (both in time and space) with no smoothing. We can see that the curve comes from the 
unstable regime, then moves to the stable one, and finally returns to the unstable regime. The whole 
simulation is carried out in the stable regime. 
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Figure 1: z{a,0) from Section[21 


3 Statement and proof of the main result 

Our main result, Theorem 13.31 below, is motivated by the above numerics (although the mechanics of 
the constructed solutions will be somewhat different). This section is devoted to its proof, which will 
rest on the following two lemmas. 


Lemma 3.1 There exists an analytic, odd, asymptotically flat initial condition z(a,0) whose analytic 
extension is on the boundary of its strip of analyticity, and also ai > 0 such that the following 
hold. We have daZi{a,0) = 2 : 2 ( 0 ;, 0) = 0 < daZ2{o',0) for a G {0,±q:i}, while daZi{a,0) > 0 for all 
other a G R (in particular, z{a, 0) is a graph of a function of x with three vertical tangents), such that 
the corresponding solution z{a,f) of (ll.4|l satisfies sgn.(t)daZi{±ai,t) > 0 and sgn.(t)daZi{0,t) < 0 
for all times t sufficiently close to 0. 

By asymptotically flat we mean 2:(a,0) « (o, 0) and daz{a,0) « (1,0) for |a| 1. Therefore, 

2 :(q:, t) will be a graph of a function of x everywhere except near a = 0 for small t > 0, and everywhere 
except near a = ± 0:1 for small t < 0. 


Proof: Since in the following we will mostly consider t = 0, let us denote z(a) = z{a,0). Following 
the arguments of [5], we find that if ao is any point with 2 :((q:o) = - 21 ( 00 ) = 2 : 2 ( 00 ) = 0, then the 
corresponding velocity v at time t = 0 satisfies 


daVi{z{ao)) 



(21 (/ 3 ) - zi{ao))z[{l3)z2{f) 

[(zflao) - zmr + Z2W)T 


We obviously have the following: 

• If daVi{z{ao)) > 0, then sgn{t)daz{a,t) > 0 for (o,t) close to (ao,0). In particular, the curve 
will be in the stable regime near ao for small t > 0. 

• If daVi{z{ao)) < 0, then sgn{t)daz{a,t) < 0 for (a,f) close to (q;o,0). In particular, the curve 
will be in the unstable regime (i.e., it will turn over) near ao for small t > 0. 
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Figure 2: Backwards evolution of z{a,t) from Section [21 with regularization: t = 0, thick line; 
t = —1.2 ■ 10“^, broken line; t = —2.4 • 10“^, dotted line; t = —3.6 • 10“^, broken and dotted line; 
t = —4.92 ■ 10“^, thin line. Note that this is a closeup of the solution near the origin. 


Our z{a) will consist of three building blocks: center and two tails. First let 


and 


zl{a) = 


z1{a) = 


(note that both are odd). 


a 

-2 < a < -1 



( -a-2 

-2 < a < -1 


-1 < a < 1 

4(«) 


a 

-1 < a < 1 

a 

1 < a < 2 



[ -a+ 2 

1 < a < 2 





' fa + T 

-7 < a < -5 





3 

-5 < a < -2 

a 

-7 < a < -1 



-5a- 7 

-2 < a < -1 


-1 < a < 1 


= < 

3a — a^ 

-1 < a < 1 

a 

1 < a < 7 



-5a-k 7 

1 < a < 2 





-3 

2 < a < 5 





^a- ^ 
^2^ 2 

5 < a < 7 

We will now consider 

curves 

z^{a) of the following form: 


zf(a) 


z{{a + R) — R 
zfia) 

z{(a — R) + R 
a 


-R-2<a<-R + 2 
-7 <a<7 
R-2<a<R+2 
otherwise 


z^{a) = < 


zl{a + R) 
Zzip) 
z\[a - R) 
0 


-i?-2<a<-i? + 2 
-7<a<7 
R-2 <a < R + 2 
otherwise 


where R > 9 will be fixed later. A sketch of the initial condition z^{a), which satisfies all hypotheses 
of the lemma (with ai = R) except regularity, appears in Figure O We will now show that for small 
t > 0, the corresponding solution of (HI) will turn over near a = 0 but not near a = ±i?. To do so, 
we will split the integrals in the formula for daVi{z^{ao)) (ao = 0,±i?) into center-center, tail-tail, 
center-tail, and tail-center contributions (note also that {z^y{aQ) > 0). We will compute the first two 
and only estimate the last two since they will be made arbitrarily small by choosing R large enough. 
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Figure 3: z^{a) from Section [31 


We start with center-center. Since is odd, the quantity we are interested in simplifies to 


Irr = 2 


zmiz^YWzm 


where 


'o 

3x^(x^ — 3) 


+ ^cc + ^cc + 4c) 


IL = -2 


ll = 2 


/ 


0 (2a;4 - 6a:2-p 9)2 
^ 7x - 5a;2 


/3 = —2 

= 


(26x2 - 70x -f 49)2 
3x 


dx 

dx 


dx 


I 


J2 (9 + x2)2 

^ 3x2 _ 2]^2; 


/441 63a; , 13a;^ \ 

V 4 2 + 4 / 


The last three integrals can be explicitly calculated: 


j 2 _ ^ f -7-flOx 


26 V 26x2 - 70x -b 49 

21=5 


It = 3 


It = 


1 


9 -b x2 


x=2 


63 

'442 


48(7 - 2x) 


338x2 - 3276x-b 11466 


x=2 


x—1 


21 = 7 


21 = 5 


dx. 


1 


3 

119' 


The first one can also be calculated explicitly, with slightly more effort. For the sake of simplicity, we 
do not present here the full symbolic integral, only an approximation of the final result: 


0.127271158. 
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Adding all the contributions, we obtain 


Ice w -0.0250882 < 0. 

Next, the tail-center term (i.e., the contribution of the tails to ( 2 :^( 0 ))) is 

/■^ {R+z\{mzi)'mzm ,0 
J-2[{R+zmY+zm^? 

The explicit expression for z* yields the easy bound 

I (i^ + 2)-3-l 

\Itc\-2 (^_ 2)4 4 

Hence, by choosing R large enough we can ensure that 

Icc + Itc < 0 . 

Let us now consider the tail-tail contribution. Because of symmetry it suffices to consider ag 
Then we have 

where (after using the explicit expression for z*) 


1 ^ ^ x{2 + x) 

“ J-2 (2a;2 +4:X + 4)2 


dx 


Iu = 




-dx 


/_i (l+a;4)2 
/3 = _ 


11 (2x2 — 4x -b 4)2 


The first and last integral can again easily be evaluated explicitly: 

X — —1 


Iu = 


T'^ — 

^tt — 


l-\- X 


4(2 + 2x + x2) 

X — \ 


4(2-2x + x2) 
sitive. We 

Finally, we bound the center-tail contribution 


x^-2 

x=2 


x—l 


whereas > 0 since its integrand is positive. We can thus conclude that Rt > \- 


where 


T _ _i_ r2 

^ct -^ct ' 


ri ^ /■" jzm - R){ziy{i3)zm 

J 2 ^ {zm-R){z\y{P)z\{P) 

J.n .2 [{R-z{{PW + zm^? 


(3.1) 
= R. 
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We can easily obtain the bounds 


1 (fl + 7)-3-3-14 

I ctl - (R-7)i 

. (2j? + 2)-3-l-4 

led- (2i?-2)4 

so again, by choosing R large enough we can ensure that 


Itt + let > 0. (3.2) 

We therefore choose R such that conditions (13.11) and (13.21) are satished. Then we let be a 

sufficiently close analytic perturbation of z^{a) which satisfies all the hypotheses of the lemma with 
ai = R (so it also has vertical tangents at a = 0, zkR). It is not difficult to see that this can be done 
so that 9aUi(2^^(ao)) has the same sign as daVi{z^{ao)) for ao = 0, ±i?. Then z{a,0) = 
will be the desired initial condition. □ 


Lemma 3.2 Let z{a,t) and w{a,t) he two analytic solutions of equation (11.41) with initial eonditions 
zo{a) and wo{a) respectively. Let d(a,f) = z{a,f) — w(a,f), and for q:,/3 G R let 


F{z,t){a,P) 


\z{a,t) - z{a- /3,t)\‘^' 


Assume F{z,0), F{w,0) G L°° and consider the energy 


E{t) = \\d{-,t)\\jji,guss = sup [ {\d{a + ic,t)\'^+ \d^d{a + ic,t)\'^) da, 

with S(t) = {a + ic \ |c| < f,(t)} the strip of analyticity of the function d{-,t). Then there exists a 
polynomial P{x,y,q,r,s) such that 

-^E{t) < P {E{t), ||2ollif4(S(o))i l|w^ollff4(S(o))i II-P'(2:,0)||loo(5(o)), ||F(ri;, 0)||ioo(s(o))) , 

where the L°° norm in the strip is defined as 

\\F{zA)\\L’^{S{t)) = sup ||i^(z,t)(a + ic,/3)||i=o(„_^). 

|c|<C(t) 

Proof: The proof appears in [3 Section 6, p. 940, eq. (44)]. One only has to write the equation that 
d{a,t) satisfies, then apply the same estimates as in the local existence theorem from [^, and hnally 
control the evolution of the norms of z and w in terms of the norms of zq and wq via their respective 
local existence theorems. □ 


We are now ready to prove our main result. 

Theorem 3.3 There exists an analytic initial curve z(a, 0) whose analytic extension is Fl'^ on the 
boundary of its strip of analyticity, and also times —T < ti < 0 < to < T and e > 0 such that the 
following hold. The corresponding solution of (113 exists for t G (—T, T) in the class of analytic 
functions of a whose analytic extensions are on the boundaries of their strips of analyticity, and 
it is a graph of a function of x for each t G (<i,to) (i.e., it is in the stable regime for these t) but not 
for t £ ifi — e,ti) U {to,to + s) (i.e., it is in the unstable regime for these t). The solution develops 
vertical tangents at times ti and to. 
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Proof: For i5 > 0, let z^{a) be a sufficiently close perturbation of 2 ^ 0 ( 0 ) with the same properties 
except that its tangents at a = 0,±i? are S away from vertical (in the sense of daz({ao) = S for 
oq = 0,±i?, while daZ^icuo) remains away from 0, uniformly in i5). It is not difficult to see that this 
can be chosen with a radius of analyticity away from 0, uniformly in i5. Then the solutions z^(a, t) of 
(11.41) corresponding to initial conditions z^{a) exist in the class of analytic functions on the interval 
(—r, T) for some (^-independent T > 0. 

Let to{6) be the time of turnover of z^{a,t) near oq = 0, and let ti{6) be the time of turnover 
near ctQ = ±R (these exist if z^{a) is close enough to 2 ^ 0 ( 0 ;)). We have to(0) = ti(0) = 0, as well 
as to{S) > 0 > ti{S) for sufficiently small 6 and lim5_>o to(<5) = lim 5 _>o ti((5) = 0, due to Lemma [321 
Choosing 6 such that —ti{S),to{6) < T yields the desired initial condition 2 ( 0 , 0) = z^{a). □ 

4 Further discussion 

We now discuss a couple of remarks related to Theorem 13.31 We plan to address these in the future. 

The first is the question of narrowing the gap between Theorems If. II and If. 21 which was our initial 
motivation. Based on our numerical simulations, we propose the following conjecture: 

Conjecture 4.1 There exists a solution f{x,t) of equation (11.51) that has ||/(•, 0)||loo = 50 and turns 
over in finite time. 

We present in Figure |T] the supporting numerical evidence, with initial condition 

2 

2 i(q;, 0) = a — 0.96sin(a), 22 ( 0 , 0) = , sin(3a). 

o 

If we reparametrize this curve as (x,/(x)), then ||/'||ioa = /'(O) = 50. 



Figure 4: /(x, 0) from Section HI Inset contains zoom (around x = 0) of f{x,t) for different times: 
t = 0, broken line; t = 5.86 • 10“^, dotted line; t = 1.26 • 10“^, thin line; t = 1.93 • 10“^, thick line. 

A possible strategy to proving this conjecture, outlined in [3], is as follows. One can consider an 
approximate solution w{a,t), satisfying (with a small error err{a,t)) 

dtw{a,t) = ^ ^ f {daw{a, t) - dpwiP, t))dp + err{a, t) (4.1) 

27r J \w[a,t) - w{l3,t)Y 
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Next, subtracting equation dSD from dEl and defining d{a,t) = z{a,t) — w{a,t), one can write 
a system of equations for d(a,t) in such a way that only d(a,t) w(a,t) and err{a,t) appear, and 
z{a,t) does not. We recall that w{a,t) is explicit since it is the numerically calculated function and 
da'Wi{a, T) < 0 for some explicit (a, T). In a similar fashion as in the local existence theorem from [^, 
a stability theorem follows, and this gives explicit bounds on the evolution of some norms of d{a^t), 
in particular it controls the L°° norm of dadi{a,t). The bounding of the constants of the stability 
theorem can be done either via traditional pencil-and-paper means, or using rigorous computer-assisted 
bounds and interval arithmetics, or a combination of both. Finally, if \\dad\{a, T)\\l<x, is small enough, 
then 9„-2 i(q:, T) = da'Wi{a,T)+dadi{a,T) <0. Note that, as opposed to Lemma l3.2l we need to have 
good enough bounds on the constants that appear in the inequality, not just to know their existence. 
This makes the problem considerably harder. 

Remark 4.2 We believe that the constant 50 is not sharp and can be improved with further numerical 
search and better estimates. 

Remark 4.3 Showing existence of a stability shift in the other direction (stable unstable —> stable) 
seems harder, even though we have numerics that exhibit that behaviour. One has to produce an initial 
condition that first turns over and then recoils back, so in contrast with Theorem \3.3\. the solution lives 
in the unstable regime during the “middle” time interval of its evolution. 

A similar strategy as outlined above could work. After finding a numerical guess for a such solution, 
one can try to show via a stability theorem and computer-assisted estimates that close to this guess 
there exists a true solution of dEl which still exhibits this phenomenon. 
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